{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2b5ce53a-e2ff-4265-ac7c-6662adca43eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "suppressMessages(library(ArchR))\n",
    "suppressMessages(library(patchwork))\n",
    "suppressMessages(library(ggplot2))\n",
    "suppressMessages(library(dplyr))\n",
    "suppressMessages(library(cowplot))\n",
    "suppressMessages(library(ggrepel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ac335ed0-5d80-4763-8710-a512b4e8ed0e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Setting default number of Parallel threads to 126.\n",
      "\n",
      "Setting default genome to Hg38.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "set.seed(42)\n",
    "addArchRThreads(threads = parallel::detectCores() - 2)\n",
    "addArchRGenome(\"hg38\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "482d46b7-4cce-4d5b-b0e5-6749571c06f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Successfully loaded ArchRProject!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "proj <- loadArchRProject(\"../data/VisiumHeart\", showLogo = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3383ed91-aa17-4310-a6f9-62ec792136c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ArchR logging to : ArchRLogs/ArchR-addGroupCoverages-3e0ae01fc58622-Date-2021-09-16_Time-11-10-03.log\n",
      "If there is an issue, please report to github with logFile!\n",
      "\n",
      "CK166 (1 of 8) : CellGroups N = 2\n",
      "\n",
      "CK167 (2 of 8) : CellGroups N = 2\n",
      "\n",
      "CK168 (3 of 8) : CellGroups N = 2\n",
      "\n",
      "CK169 (4 of 8) : CellGroups N = 2\n",
      "\n",
      "CK170 (5 of 8) : CellGroups N = 2\n",
      "\n",
      "CK171 (6 of 8) : CellGroups N = 2\n",
      "\n",
      "CK173 (7 of 8) : CellGroups N = 2\n",
      "\n",
      "CK174 (8 of 8) : CellGroups N = 2\n",
      "\n",
      "2021-09-16 11:10:08 : Creating Coverage Files!, 0.09 mins elapsed.\n",
      "\n",
      "2021-09-16 11:10:08 : Batch Execution w/ safelapply!, 0.09 mins elapsed.\n",
      "\n",
      "2021-09-16 11:10:58 : Adding Kmer Bias to Coverage Files!, 0.924 mins elapsed.\n",
      "\n",
      "Completed Kmer Bias Calculation\n",
      "\n",
      "Adding Kmer Bias (1 of 16)\n",
      "\n",
      "Adding Kmer Bias (2 of 16)\n",
      "\n",
      "Adding Kmer Bias (3 of 16)\n",
      "\n",
      "Adding Kmer Bias (4 of 16)\n",
      "\n",
      "Adding Kmer Bias (5 of 16)\n",
      "\n",
      "Adding Kmer Bias (6 of 16)\n",
      "\n",
      "Adding Kmer Bias (7 of 16)\n",
      "\n",
      "Adding Kmer Bias (8 of 16)\n",
      "\n",
      "Adding Kmer Bias (9 of 16)\n",
      "\n",
      "Adding Kmer Bias (10 of 16)\n",
      "\n",
      "Adding Kmer Bias (11 of 16)\n",
      "\n",
      "Adding Kmer Bias (12 of 16)\n",
      "\n",
      "Adding Kmer Bias (13 of 16)\n",
      "\n",
      "Adding Kmer Bias (14 of 16)\n",
      "\n",
      "Adding Kmer Bias (15 of 16)\n",
      "\n",
      "Adding Kmer Bias (16 of 16)\n",
      "\n",
      "2021-09-16 11:11:39 : Finished Creation of Coverage Files!, 1.608 mins elapsed.\n",
      "\n",
      "ArchR logging successful to : ArchRLogs/ArchR-addGroupCoverages-3e0ae01fc58622-Date-2021-09-16_Time-11-10-03.log\n",
      "\n"
     ]
    }
   ],
   "source": [
    "proj <- addGroupCoverages(ArchRProj = proj, groupBy = \"Sample\",\n",
    "                          maxCells = 5000,\n",
    "                          minCells = 5000,\n",
    "                          force = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cd170871-30ae-4e23-bedc-f19b38b099ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Searching For MACS2..\n",
      "\n",
      "Found with $path!\n",
      "\n",
      "ArchR logging to : ArchRLogs/ArchR-addReproduciblePeakSet-3e0ae011b57349-Date-2021-09-16_Time-12-07-05.log\n",
      "If there is an issue, please report to github with logFile!\n",
      "\n",
      "Calling Peaks with Macs2\n",
      "\n",
      "2021-09-16 12:07:05 : Peak Calling Parameters!, 0.003 mins elapsed.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Group nCells nCellsUsed nReplicates nMin nMax maxPeaks\n",
      "CK166 CK166   3402        540           2   40  500   150000\n",
      "CK167 CK167   2575        540           2   40  500   150000\n",
      "CK168 CK168   2156        540           2   40  500   150000\n",
      "CK169 CK169    856        540           2   40  500   150000\n",
      "CK170 CK170   1651        540           2   40  500   150000\n",
      "CK171 CK171   5033        540           2   40  500   150000\n",
      "CK173 CK173   1743        540           2   40  500   150000\n",
      "CK174 CK174    284        284           2   40  244   142000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-09-16 12:07:06 : Batching Peak Calls!, 0.003 mins elapsed.\n",
      "\n",
      "2021-09-16 12:07:06 : Batch Execution w/ safelapply!, 0 mins elapsed.\n",
      "\n",
      "2021-09-16 12:10:21 : Identifying Reproducible Peaks!, 3.258 mins elapsed.\n",
      "\n",
      "2021-09-16 12:10:33 : Creating Union Peak Set!, 3.464 mins elapsed.\n",
      "\n",
      "Converged after 6 iterations!\n",
      "\n",
      "Plotting Ggplot!\n",
      "\n",
      "2021-09-16 12:10:41 : Finished Creating Union Peak Set (223167)!, 3.595 mins elapsed.\n",
      "\n",
      "ArchR logging to : ArchRLogs/ArchR-addPeakMatrix-3e0ae010664d98-Date-2021-09-16_Time-12-10-41.log\n",
      "If there is an issue, please report to github with logFile!\n",
      "\n",
      "2021-09-16 12:10:41 : Batch Execution w/ safelapply!, 0 mins elapsed.\n",
      "\n",
      "ArchR logging successful to : ArchRLogs/ArchR-addPeakMatrix-3e0ae010664d98-Date-2021-09-16_Time-12-10-41.log\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## peak calling\n",
    "pathToMacs2 <- findMacs2()\n",
    "\n",
    "proj <- addReproduciblePeakSet(\n",
    "    ArchRProj = proj, \n",
    "    groupBy = \"Sample\", \n",
    "    cutOff = 0.05,\n",
    "    pathToMacs2 = pathToMacs2\n",
    ")\n",
    "\n",
    "proj <- addPeakMatrix(proj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a24a4082-d7ee-4043-8118-b2df83c41fe4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ArchR logging to : ArchRLogs/ArchR-getMatrixFromProject-3e0ae0294c2410-Date-2021-09-16_Time-12-14-07.log\n",
      "If there is an issue, please report to github with logFile!\n",
      "\n",
      "2021-09-16 12:14:40 : Organizing colData, 0.55 mins elapsed.\n",
      "\n",
      "2021-09-16 12:14:40 : Organizing rowData, 0.552 mins elapsed.\n",
      "\n",
      "2021-09-16 12:14:40 : Organizing rowRanges, 0.552 mins elapsed.\n",
      "\n",
      "2021-09-16 12:14:41 : Organizing Assays (1 of 1), 0.552 mins elapsed.\n",
      "\n",
      "2021-09-16 12:14:55 : Constructing SummarizedExperiment, 0.791 mins elapsed.\n",
      "\n",
      "2021-09-16 12:14:56 : Finished Matrix Creation, 0.804 mins elapsed.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## save peak matrix\n",
    "peakMatrix <- getMatrixFromProject(proj,\n",
    "                                   useMatrix = \"PeakMatrix\")\n",
    "\n",
    "counts <- peakMatrix@assays@data$PeakMatrix\n",
    "df_rangers <- as.data.frame(peakMatrix@rowRanges@ranges)\n",
    "\n",
    "rownames(counts) <- paste(peakMatrix@rowRanges@seqnames,\n",
    "                          df_rangers$start,\n",
    "                          df_rangers$end,\n",
    "                          sep = \"_\") \n",
    "\n",
    "saveRDS(counts, file = \"../data/VisiumHeart/PeakMatrix.Rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6e977d1c-a463-41a8-92ba-42462205a36c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ArchR logging to : ArchRLogs/ArchR-getMatrixFromProject-3e0ae04cfea0c5-Date-2021-09-16_Time-12-17-32.log\n",
      "If there is an issue, please report to github with logFile!\n",
      "\n",
      "2021-09-16 12:17:58 : Organizing colData, 0.434 mins elapsed.\n",
      "\n",
      "2021-09-16 12:17:58 : Organizing rowData, 0.436 mins elapsed.\n",
      "\n",
      "2021-09-16 12:17:58 : Organizing rowRanges, 0.436 mins elapsed.\n",
      "\n",
      "2021-09-16 12:17:58 : Organizing Assays (1 of 1), 0.436 mins elapsed.\n",
      "\n",
      "2021-09-16 12:18:14 : Constructing SummarizedExperiment, 0.709 mins elapsed.\n",
      "\n",
      "2021-09-16 12:18:15 : Finished Matrix Creation, 0.724 mins elapsed.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## save gene matrix\n",
    "atac <- getMatrixFromProject(ArchRProj = proj,\n",
    "                             useMatrix = \"GeneScoreMatrix\")\n",
    "\n",
    "counts <- atac@assays@data$GeneScoreMatrix\n",
    "rownames(counts) <- atac@elementMetadata$name\n",
    "\n",
    "saveRDS(counts, file = \"../data/VisiumHeart/GeneScoreMatrix.Rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72819eff-55c0-4751-a3a4-244b9b830f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "saveArchRProject(ArchRProj = proj, \n",
    "                 load = FALSE)\n",
    "sessionInfo()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
